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Abstract 

The Banks-Casher relation links the spontaneous breaking of chiral symmetry in QCD to 
the presence of a non-zero density of quark modes at the low end of the spectrum of the 
Dirac operator. Spectral observables like the number of modes in a given energy interval 
are renormalizable and can therefore be computed using the Wilson formulation of lattice 
QCD even though the latter violates chiral symmetry at energies on the order of the inverse 
lattice spacing. Using numerical simulations, we find (in two-flavour QCD) that the low 
quark modes do condense in the expected way. In particular, the chiral condensate can be 
accurately calculated simply by counting the low modes on large lattices. Other spectral 
observables can be considered as well and have a potentially wide range of uses. 



1. Introduction 

So far all results obtained in numerical lattice QCD are consistent with the expec- 
tation that chiral symmetry is spontaneously broken in the way presumed by chiral 
perturbation theory. Little is known, however, about the dynamical processes that 
cause the symmetry to break. An intriguing remark, made long ago by Banks and 
Casher [1], is that the effect is tied to a condensation of the low modes of the Dirac 
operator. Studies of the low modes may therefore provide important clues on the 
symmetry-breaking mechanism. 

In the Wilson formulation of lattice QCD [2] and its improved versions [3,4], chiral 
symmetry is violated explicitly by terms proportional to the first or second power of 
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the lattice spacing. The Banks-Casher relation consequently cannot be expected to 
hold exactly and the detailed properties of the low quark modes could be significantly 
different from those in the continuum theory. On the other hand, as long as only 
renormalizable quantities are considered, their values in the continuum limit must 
in principle be computable using the Wilson theory. 

The spectral density of the (hermitian) Dirac operator, and thus the average num- 
ber of quark modes in a given range of eigenvalues, are known to be renormalizable 
[5]. In the present paper, we first give a second proof of this important fact (sect. 3). 
We then discuss the chiral perturbation expansion of the mode numbers and show, in 
sect. 5, that their calculation in lattice QCD requires only a modest computational 
effort. Taken together, these results allow the chiral condensate to be computed in 
the Wilson theory in a straightforward manner (sect. 6). Spectral projectors however 
have a wider range of applicability and provide interesting opportunities to explore 
the chiral regime of QCD, some of which are briefly mentioned in sect. 7. 



2. Preliminaries 

For simplicity we focus on QCD with a doublet of mass-degenerate quarks, but the 
theoretical discussion is more generally valid and extends to the case of real-world 
QCD. The quarks will be referred to as the up and down quarks, the associated 
Goldstone bosons as the pions and the SU(2) flavour symmetry as the isospin sym- 
metry. We consider both the continuum and the Wilson lattice theory in order to 
make it clear in which way the mode number computed on the lattice is related to 
the one defined in the continuum theory. 

2.1 Spectral density and mode number in the continuum theory 

In a space-time box of volume V with periodic or antiperiodic boundary conditions, 

the euclidean massless Dirac operator D in presence of a given gauge field has purely 
imaginary eigenvalues iAi, 1X2, . . which may be ordered so that those with the 
lower magnitude come first. The associated average spectral density is given by 



where the bracket (. . .) denotes the QCD expectation value and m the current-quark 




k=i 



(2.1) 
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mass. Note that the isospin degeneracy is not included in the mode counting, i.e. the 
Dirac operator is diagonalizcd in the subspace of, say, the up-quark fields. 
The Banks-Casher relation [1] 

lim lim lim p(A, m) = — 

provides a link between the chiral condensate 
S = — lim lim (uu) 

m— >0 V—>oo 

(where u is the up-quark field) and the spectral density. In particular, if chiral sym- 
metry is spontaneously broken by a non-zero value of the condensate, the density of 
the quark modes in infinite volume does not vanish at the origin. A non-zero density 
conversely implies that the symmetry is broken, i.e. the Banks-Casher relation can 
be read in either direction. 

Instead of the spectral density, the average number i^(M, m) of eigenmodes of the 
massive hermitian operator D^D + 171^ with eigenvalues a < turns out to be a 
more convenient quantity to consider. Evidently, since 

/.A 

u{M,m) = V dXp{X,m), A= VM^-m^, (2.4) 

J-A 

the mode number ultimately carries the same information as the spectral density. 
2.2 0(a) -improved lattice QCD 

The lattice theory is set up as usual on a hyper-cubic lattice with spacing a, time-like 
extent T and spatial size L. Periodic boundary conditions are imposed on all fields 
and in all directions, the only exception being the quark fields which are taken to 
be antiperiodic in time. 

As already mentioned, we focus on the Wilson theory in this paper. The details 
are not very relevant, but for definiteness we choose the Wilson plaquette action for 
the gauge field [2] and the standard expression 

Sp = a'^^ {u{x)Dmu{x) + d{x)Dmd{x)} (2.5) 

X 

for the quark action, in which Dm denotes the massive, 0(a)-improved lattice Dirac 
operator [3,4]. Apart from the bare coupling gg and the bare mass mo, the only free 
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(2.2) 



(2.3) 



parameter in the lattice action is the improvement coefficient Cgwj which we choose 
so as to cancel the 0(a) lattice effects in on-shell quantities [6]. 

In this theory, the renormalized coupling and quark mass are related to the bare 
parameters through [4] 

= Zg{l + hgam^)gl, (2.6) 
= Zmi'i- + bmamqjniq, mq = mo - rric, (2.7) 

where mc{go) denotes the critical bare mass and bg{go) and bm{go) are further 0(a)- 
improvement coefficients. The renormalization constants Zg and depend on the 
normalization conditions and are functions of the bare coupling and a normalization 
scale given in units of the lattice spacing. 

Composite fields like the isospin axial current and the isospin pseudo-scalar and 
scalar densities are renormalized similarly by factors of the form Zx{l + bxaniq) 
where X = A,P, S. The normalization conditions will be assumed to be such that 
the renormalized correlation functions satisfy the non-singlet chiral Ward identities 
up to terms of order a^. In particular, 

ZA{l + bAamq) 2n /o 

"^"= Zp(l + 6pamj "^ + Q^" ^' ^^-'^ 

where m is the bare current-quark mass that appears in the PCAC relation [4]. 

On the lattice, we shall be mostly interested in the average number niq) 
of eigenmodes of D„j D„i with cigcnvahics a < Af^. This definition of the mode 
number formally coincides with the one given in subsect. 2.1, but it would evidently 
be premature to conclude that the values calculated on the lattice are simply related 
to the mode number defined in the continuum theory. 



3. Renormalization of the mode number 



The proof of the renormalizability of the mode number given in this section partly 
follows the lines of ref. [5] , but avoids some of the rather technical assumptions that 
had to be made there. An important new element of the proof is the use of twisted- 
mass valence quarks and the associated density-chain correlation functions, which 
have other applications as well (see sect. 7). 
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3.1 Spectral sums and density chains 

We consider the lattice theory and introduce the spectral sums 

<7fc(M,mq) = (TV{(i^Jl)^ + /x2)-'=}), (3.1) 

where A; > 3 will be assumed for reasons to become clear below. The spectral sums 
are related to the mode number ^{M, niq) through the integral transform 

f°° 2kM 
akili,m^)= dMv{M,m^) (3.2) 
Jo [M'^ + ji^) 

which can be shown to be invertible for every fixed k. The renormalization properties 
of z/(M, mq) can therefore be inferred from those of, say, a^{ii,mq). 

The inverse of the operator D^D^ + /x^ coincides with the square of the quark 
propagator in twisted-mass lattice QCD [7]. We are thus led to add a set of isospin 
doublets ipi, I = 1, . . . , 2k, of valence-quark fields to the theory, with action 

2k 

^ ^ Mx) {Dm + Wsr^) M^) (3.3) 

X 1=1 

(the isospin indices are suppressed in this formula and is the third isospin Pauli 
matrix). Evidently, in order to cancel the valence-quark determinant, a correspond- 
ing multiplet of pseudo-fermion fields must be added as well. The spectral sums 
(3.1) can then be represented by density-chain observables like 

(j3(//,m) = -a^^ X] 

Xi,...,X6 

(P+(a;i)P2-3(x2)P3+4(a;3)P45(^4)P5'6(^5)^'6"l(^6)>, (3.4) 

where = V'iTsT^Vj the charged pseudo-scalar densities of the valence quarks 
(see fig. 1). 

3.2 Renormalization of the spectral sums 

With respect to the case of twisted-mass QCD discussed by Prezzotti et al. [7,8], the 
0(a)-improvement and renormalization of the partially quenched theory considered 
here tends to be somewhat simpler. In particular, we may choose a scheme which is 
independent of the twisted mass parameter and which coincides with the commonly 
used conventions in the sea-quark sector of the theory. 
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Fig. 1. The flavour labels of the pseudo-scalar densities in eq. (3.4) are such that 
the contraction of the quark fields yields a closed quark loop with six edges. Each edge 
represents a propagator {Dm ± i/Lt75)~^ and each vertex contributes a factor 75. The 
ordered product of these factors, summed over the positions xi, . . . ,X6 of the fields, 
coincides with the trace (3.1). 

At n = 0, the Wilson theory preserves the lattice symmetries, charge conjugation, 
the gauge symmetry and all (vector) flavour symmetries, including the ones that 
mix the sea with the valence quarks. Ultraviolet-divergent terms other than those 
cancelled by the usual parameter and field renormalizations are excluded by these 
symmetries. When the twisted mass /j, is switched on, some of the symmetries are 
broken and further ultraviolet-divergent terms can arise. Power counting then shows 
that a multiplicative renormalization, 



plus the renormalizations required at )U = are sufficient to renormalize the partially 
quenched theory. Moreover, the correction proportional to arriq included in eq. (3.5) 

is all what needs to be added for on-shell 0(a)-improvemcnt at /U 7^ [8]. 

Considering eq. (3.4), these remarks suggest that the renormalization of cr^lp, mq) 
is achieved by multiplication with the sixth power of the renormalization factor Zp 
of the pseudo-scalar densities and by renormalizing the parameters of the theory. 
The only worry one may have at this point is that the summations in eq. (3.4) 
over the coordinates xi, . . . ,xq diverge in the continuum limit. However, as already 
pointed out in refs. [9,10,5], the short-distance singularities of density-chain correla- 
tion functions are integrable, and give rise to O(amq) corrections only, if there are 
six or more densities. 

For any k > S, the renormalized 0(a)-improved spectral sums are thus given by 



Mr = Zij,{l + b^amq)n, 



(3.5) 



1 -|- bparriq 



} 



2k 



Zp 



1 -|- bppaniq 



(Tfe(|U,mq), 



(3.6) 
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where it is understood that the bare masses are expressed through the renormahzed 
ones. The factors 1 + bppaniq in eq. (3.6) are required for the cancellation of the 
O(amq) terms alluded to above which derive from the short-distance singularities of 
the density-chain correlation functions [5]. 

3.3 Renormalized mode number 

If the twistcd-mass term is considered to be a perturbation of the theory at /x = 0, 
one quickly notices that 

= (3-7) 

is a possible (and natural) choice of the renormalization factor Z^. 
Another simplification derives from the identity 

d 

— cjfc(/i,mq) = -2A;Aicrfc+i(/x,mq). (3.8) 

When the renormalized spectral sums are similarly differentiated with respect to the 
renormalized twisted mass /Xr, the expressions one obtains must be 0(a)-improved. 
As it turns out, this is the case if and only if 

b^, + bp- bpp = 0. (3.9) 

The renormalization factor in eq. (3.6) thus becomes 

1 + bpam^ ^ 1 
1 -I- bppam^ Zn{l + b^j^am^ 

up to terms of order a^rriq. 

Returning to the integral representation (3.2), we now note that the renormaliza- 
tion factor {Zfj,{l + bfj,amq)}~^'^ needed to renormalize the spectral sum on the left 
of the equation is cancelled on the right if we substitute 

Mr = Z^(l-F6^amq)M (3.11) 

and renormalize /x. We are thus led to conclude that 

MMr, mn) = u{M, mq) (3.12) 

is a renormalized and 0(a)-improved quantity. In other words, the mode number is 
a renormalization-group invariant. 
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3-4 Universality 

The steps taken in this section can be repeated using other regularizations of QCD 
as long as these preserve the same (or more) symmetries as the Wilson theory. Di- 
mensional regularization with the 't Hooft-Veltman prescription for 75, for example, 
has all the required properties, although in this case one is limited to weak-coupling 
perturbation theory. 

Independently of the regularization, the renormalized mode number will be the 
same if the same normalization conditions are used. In particular, a definite conven- 
tion such as the MS scheme must be adopted for the normalization of the pseudo- 
scalar densities. The normalization of the sea-quark mass rn-R is then determined by 
the PCAC relation, while the one of /Xr is fixed by requiring the identity 

d 

o-fc,R(A«R, rnji) = -2A;/XR(Tfc+i,R(/XR, m^) (3.13) 



to hold after removal of the regularization. At this point, the renormalized spectral 
sums are uniquely determined and so is the renormalized mode number, since the 
integral transform 

/"°° 2kM 
o-fc,R(MR, "^r) = / dMRi/R(MR,mR) ^ (3.14) 

Jo {Ml + 

is free of normalization ambiguities. 



4. Chiral expansion of the mode number 



In the continuum theory and for small masses, the mode number can be calculated 
analytically in chiral perturbation theory. Although all results quoted below are for 
the renormalized mode number, we omit the subscript "R" in this section in order 
to simplify the notation. 

4-1 Chiral perturbation theory 

At present the chiral expansion of the spectral density p{X, m) is known to next-to- 
leading order of chiral perturbation theory. The first computation to this order was 
performed by Smilga and Stern [11] in the massless theory in infinite volume. Later 
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Osborn et al. [12] and Damgaard et al. [13] performed a more complete and system- 
atic computation based on partially quenched chiral perturbation theory [14,15]. 
The starting point in the paper of Osborn et al. is the formula 

p(A, m) = lim {Svai(e + iX) + Svai(e - iX)} , (4.1) 

which relates the spectral density to the expectation value — Svai(?TZvai) of the scalar 
density of an added valence quark of mass mvai- With a doublet of sea quarks, the 
relevant graded flavour symmetry group is then SU(3|1) and the chiral expansion of 
Svai("^vai) is derived from the associated chiral effective theory (see appendix A). 

4-2 Large-volume regime 

In infinite volume, chiral perturbation theory yields an expansion of p{X, m) essen- 
tially in powers of A and m. The leading-order term is given by the Banks-Casher 
formula and the "eflFective chiral condensate" , defined through 

therefore coincides with S in the chiral limit. 

At next-to-leading order, the chiral expansion reads 



Seff 



mS f , AS „- , -, / "T-^X 



m, A A m\ 

+ — arctan 1 arctan — > -|- . . . (4.3) 

A mm A J 

The constants F and in this expression are, respectively, the pion decay constant 
in the chiral limit and an SU(3jl) low-energy effective coupling renormalized at scale 
Ji (appendix A). Following the tradition [16], Jl may be set to the physical charged- 
pion mass, but since only the scale-invariant sum of the first two terms in the curly 
bracket matters, this choice is not compulsory. 

A remarkable feature of eq. (4.3) is that the one-loop correction vanishes, for any 
value of A, when the quark mass goes to zero. Smilga and Stern [11] already noted 
the absence of terms proportional to A and showed that this was a special property 
of the two-fiavour theory. The chiral corrections to Sefif/S consequently tend to be 
quite small (see fig. 2 for illustration). 
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Fig. 2. Quark-mass dependence of Eeflp/S at fixed A according to next-to-leading 
order of chiral perturbation theory. The low-energy constants have been set to S = 
(250MeV)3, F = 90 MeV, fi = 140 MeV and le = 3 in this plot. 

4-3 Finite-volume effects 

In the present context, the kinematical situation of interest is the so-called p-regime 
of QCD, where T > L, FL > 1 and mT,V ^ F^LF' . Chiral perturbation theory is 
easily extended to this regime and can be used to estimate the effects of the finite 
volume [17]. 

In the case of Sefr, the calculation shows that the dependence on the volume sets 
in at one-loop order and that the infinite-volume limit is reached at an exponential 
rate according to 

2eff-Seff|^_oce-*^-^ Ml = ^. (4.4) 

Note that Ma coincides with the leading-order expression for the mass of a pseudo- 
scalar meson made of two valence quarks of mass A. Since A is normally taken to be 
significantly larger than the sea-quark mass, the finite-size effects (4.4) tend to be 
smaller than those expected for the pion mass M^^, for example, which decrease like 
^-M„L^ In particular, if the parameter values previously used in fig. 2 are inserted, 
and if L > 2 fm is assumed, Ses is estimated to deviate from its infinite- volume 
value by a fraction of percent at most. 
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5. Counting the low modes in lattice QCD 



In presence of a given gauge field, the number of eigenmodes of Dm with eigenval- 
ues a < can be determined straightforwardly by calculating the eigenvalues and 
their multiplicities numerically. The effort required for such computations however 
grows proportionally to the second or perhaps even a higher power of the space-time 
volume V. In this section, we show that the modes can be counted more efficiently 
using spectral projectors. 

5. 1 Stochastic representation of the mode number 

Let Pm be the orthogonal projector to the subspace of quark fields spanned by the 
eigenmodes of Dj,}D^ with eigenvalues a < M^. An alternative representation of 
the mode number 



z/(M,mq) = (TV{Pm}) (5.1) 
is then given by 

TV 

N 



1 ^ 

u{M,m^) = {On), ^'jv = T7 (%,Pm%) , (5.2) 



fc=i 



where we have added a set of pseudo-fermion fields, iji, . . . ,r)N, to the theory with 
action 



N 



Sn = ^{Vk,Vk)- (5.3) 
fc=i 

In the course of a numerical simulation, these fields are generated randomly, for each 
gauge-field configuration, and the mode number is estimated in the usual way by 
averaging the observable On over the generated ensemble of fields. 
The variance of On, 

{{On - {On))^) = <(Tr{PM} - (TV{Pm}))'> + ^^^(M, mj, (5.4) 

is larger than the one of TrlP^}) but the difference can be reduced by increasing the 
number N of pseudo-fermion fields. More important may be the fact that the mode 
number is an extensive quantity, while the variance of TrjP^} does not appear to 
grow with the volume V of the lattice at the values of M of interest [18]. At fixed N 
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and for a given statistics, the relative statistical error of the calculated mode number 
is therefore expected to decrease like V~^l'^. 

5.2 Rational approximation 

The projector Fm can be approximated fairly easily by rational functions of DjJ Dm- 
There are different ways to proceed and the choices made in the following may not 
be the best ones, but the proposed method is quite efficient and numerically safe. 
Let P{y) be the minmax polynomial of degree n which minimizes the deviation 

<5= max 11- (5.5) 

e<y<l 

The numerical computation of this polynomial for specified values of n and e > is 
a standard task in approximation theory (see ref. [19], for example). In the range 
— 1 < X < 1, the function 

h{x) = ^ {1 - xP{x^)} (5.6) 

then provides an approximation to the step function 6(—x). By construction, the 
approximation error is at most ^5 if |x| > ^/e and numerical inspection moreover 
shows that h{x) decreases monotonically in the transition region \x\ < -y/e. 
An approximation to the projector Pjv^ is now given by 

Fm ^ h(X)\ X = 1 ; i , (5.7) 

where ~ M is an adjustable mass parameter. The quality of the approximation 
is determined by the values of n, e and the ratio M/M*. In practice, the degree n 
of the minmax polynomial should be reasonably small and the deviation 

A = (TV{Pm - h{X)^]) (5.8) 

must be much smaller than the statistical errors of the calculated mode numbers. 

The estimation of A and the choice of M/M* are discussed in appendix B. Here 
we only note that the computation of 

{'njMv) ivM^fn) = ll^i Wr/f (5.9) 

requires the application of the square of /i(X) to the pseudo-fermion field rj and not 
of its fourth power. 
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5.3 Numerical implementation 

The minmeix polynomial P{y) and therefore the operator h{%) can be expanded in a 

series of Chebyshcv polynomials with rapidly decreasing coefficients [19]. Chebyshev 
scries of this kind can be safely evaluated using the Clcnshaw recursion [20] . 

The computation of h{%)r] for a given source field r] then requires the operator 
X to be applied 2n + 1 times. Each application essentially amounts to solving the 
linear system 

{DjDm + Ml)i: = ri (5.10) 

using one's favourite iterative algorithm. This system is normally significantly better 
conditioned than the lattice Dirac equation Dmip = if]. Moreover, it is our experi- 
ence that a fairly loose stopping criterion can be chosen without compromising the 
correctness of the simulation results. 

We finally remark that the computational effort required for the calculation of the 
mode number along the lines explained here scales like V or at most V ln(F) as the 
lattice is increased. 



6. Computation of the chiral condensate 

The simulations discussed in this section have a limited scope, but the results clearly 
show that the low modes of the Dirac operator condense and that the mode number 
can be accurately computed using the stochastic method described in the previous 

section. 

We have considered two lattices in these studies, with spacing a ~ 0.08 fm, spatial 
sizes L ~ 1.9 fm and 2.5 fm, respectively, and time-like extents T = 2L. The exact 
parameter values and further technical details are given in appendix C. All values 
quoted for the renormalized mass parameters, the mass-dependent condensate Sr 
defined in subsect. 6.3 and the condensate S refer to the MS scheme at 2 GeV. 

6. 1 Qualitative behaviour of the mode number 

The data plotted in fig. 3 show that the mode number is, in the case considered, a 
nearly linear function of Mr from above the threshold region at Mr 2± ttzr up to 
at least 110 McV. This behaviour is qualitatively in line with chiral perturbation 
theory, but the fact that the linear regime extends to such large values of Mr is 
rather striking and could not be foreseen. 
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Fig. 3. Dependence of the renormalized mode number on Mr at mn ~ 26 MeV and 
L ~ 2.5 fm. The curve shown is based on a representative ensemble of 71 gauge-field 
configurations and required the lowest 80 eigenvalues of Drrl Dm to be calculated for 
each of these fields. Statistical errors are slightly larger than the jitter of the curve. 

At the very low end of the spectrum, the curve shown in fig. 3 however clearly 
deviates from its expected form in the continuum theory (shaded area in fig. 3) [5] . A 
plausible explanation of the observed deviation is that chiral symmetry is not exactly 
preserved in the Wilson theory and that the fine structure of the spectrum of the 
Dirac operator near the threshold at Mr = ttir is consequently not protected from 
perturbing lattice effects [21]. The deviation must in any case be a lattice artefact, 
since the renormalized mode number is bound to converge to its continuum value 
as the lattice spacing is decreased (cf. sect. 3). 

In the following, we focus on the linear regime in fig. 3, where the mode number 
is not expected to be particularly sensitive to discretisation errors. Moreover, since 
the effort required for the numerical calculation of the low eigenvalues of Dm is 
not small, the mode number was normally computed using the method described in 
sect. 5 and we shall, from now on, only discuss results obtained in this way. 

6.2 Volume- dependence of the mode number 

In the large- volume regime of the theory, z^(M, mq)/V is expected to be independent 
of the lattice size up to exponentially small corrections (cf. sect. 4). The lattices we 
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Fig. 4. Simulation results for the renormalized mode number at fixed L ~ 2.5 fm 

(plot on the left). The linear extrapolation to the chiral limit (open square) of Sr at 
Mr = 95 MeV is shown on the right. All errors in these plots are statistical only. 



have simulated are such that we can immediately check whether these corrections 
are significant at the level of the statistical errors. 
To this end, we form the ratios 

_ t/(M,mq)D3 /32y _ t/(M,mq)D, /32y 

""^'^ uiM,m^h, [24) ' "^''^ v{M,m^)^^ ' ^^"'^ 

where the subscripts etc. refer to the run label quoted in table 2 (appendix C). 
Both ratios turn out to be practically equal to 1. More precisely, r3^4 differs from 1 
by —0.6 to —2.0 standard deviations and rs^s by +0.7 to +1.5 standard deviations 
as M varies over the values listed in table 2. There are thus no indications for 
significant finite- volume effects on these lattices. 

6. 3 Calculation of S 

The values of the renormalized mode number which we calculated on the larger of 
the two lattices considered are plotted in fig. 4 (left graph). At fixed quark mass, the 
mode number is, to a very good approximation, a linear function of Mr in the range 
shown in the figure. In particular, the slope of the data can easily be determined by 
quadratic interpolation (lines in the left graph). 
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Table 1. Simulation results for Sr at Mr = 95 MeV 



Run 



TTiR [MeV] 



[MeV] 



45.8(3)(11) 

26.5(2)(6) 

12.8(2)(3) 



310(2)(4) 
295(2) (4) 
286(2) (4) 



We are thus led to introduce the mass-dependent condensate 




twrN ^ d 



I/R(MR,mR), 



(6.2) 



where the prefactor is chosen such that Sr coincides with the chiral condensate S to 
leading order of chiral perturbation theory. In table 1 we list the calculated values 
of Sr at Mr = 95 MeV (a point in the middle of the available range of masses). 
The first errors quoted in the table are the statistical ones, while the second errors 
are those inherited from the product of the lattice spacing and the renormalization 
factors needed to convert from lattice to physical normalizations (appendix C). 

The next-to-leading order formula (4.3) suggests that Sr = S up to higher-order 
corrections and terms vanishing proportionally to mR in the chiral limit. Note that 
there are no terms proportional to mRlnrriR at this order of the chiral expansion. 
The data for Er at Mr = 95 MeV actually fall on a straight line (right graph in 
fig. 4) and the extrapolation to mR = then yields the estimate 



for the chiral condensate. Higher-order corrections were neglected here, but appear 
to be small as the results vary only little (within roughly the third error given above) 
when the chiral limit is taken at other values of Mr. 

It goes without saying, however, that this procedure and the quoted result for the 
condensate will have to be confirmed by more extensive calculations. Meanwhile we 
note that the estimate (6.3) is in the range of values obtained in two- and three- 
flavour QCD from chiral fits of the quark-mass dependence of the pion mass [22-24] 
and from studies of the so-called e-regime of QCD [25-29]. 



SV3 ^ 276(3) (4) (5) MeV 



(6.3) 
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7. Further uses of spectral observables 



Spectral observables like the mode number provide interesting probes of low-energy 
QCD. In this section we wish to show that the computation of the chiral condensate 
is only one of the possible applications of these observables. 

7.1 Scaling to the continuum limit 

Extrapolations to the continuum limit require simulations of a series of lattices with 
decreasing lattice spacings. Since only the bare coupling and bare quark mass can 
be prescribed, the ratios of the spacings of the simulated lattices are not known a 
priori and need to be calculated. Evidently, it is very important to obtain the ratios 
with small statistical and systematic errors. 

A set of 0(a)-improved renormalized quantities, which may conceivably be used 
to match the lattices, is f 

{m,,MrG,,r,^}, (7.1) 

where M^^ and G7r,R are, respectively, the pion mass and the renormalized vacuum- 
to-pion matrix element of the isospin pseudo-scalar density. All these quantities are 
renormalization-group invariants. In particular, the dimensionless combinations 



V 



1/2 

C, = Ml[j-] , (7.2) 



3 



C2 = (MrG,,r)'' j , (7.3) 

are well-defined and directly accessible functions of go,amQ and aM. 

Since Ci and C2 are roughly linearly rising with arriQ and aM, respectively, it is 
possible to match the mass parameters on any given pair of lattices by requiring 
Ci and C2 to assume the same (sensibly chosen) values. After that the ratio of the 
lattice spacings is obtained through 



0-2 \n1V2 



f The list of observables given here only serves to illustrate the general ideas. In particular, the 
combination Mj^Er may be used in place of ui^^/V . 
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where ni , n2 denote the numbers of points of the two lattices and vi , 1^2 the mode 
numbers at the matched values of the mass parameters (we implicitly assumed here 
that finite- volume effects can be neglected or that the volumes are the same). 

An important technical advantage of this procedure is that all quantities involved 
are easily obtained with small errors. In particular, the statistical precision that can 
be attained in practice is not expected to change dramatically as the lattice spacing 
is decreased or if larger lattices are considered. 

7. 2 Computation of renormalization constants 

Density-chain correlation functions like the ones discussed in sect. 3 satisfy various 
chiral Ward identities in the continuum limit. We may, for example, start from the 
"twisted spectral sums" 

(7fe,;(/x,mq) = {Tr{j5{DjD^ + nY''j^{DjDm + fiY'})^ (7-5) 
which can be represented through density-chain correlation functions of the form 

'7l,2(/^,"Zq) = -a^^ ^ 

(5+(a;i)P23(x2).S3+(x3)P4"5(^4)P5'6(^5)P6-i(x6)>. (7.6) 

In the continuum limit and if A; -|- Z > 3, chiral symmetry (or simply the fact that 
75 commutes with Drr^ Dm in the continuum theory) then implies that the properly 
renormalized twisted spectral sum (Jk,i,K coincides with ak+i,R- 

On the lattice one should keep track of the 0(a) corrections, but following the 
lines of ref. [5] , it is then straightforward to show that 

^ = (l + 26flamq)^ + 0(a2), (7.7) 

bR = bs-bp + 2{bpp-bps), (7.8) 

where the improvement coefficient bji is known to one-loop order of perturbation 
theory and appears to be small (appendix C). 

Equation (7.7) is actually a special case of a more general Ward identity, where 
the inverse powers of D^, + the definition of the spectral sums are replaced 
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by any sufficiently rapidly decaying functions of D^D^. In particular, 

= (1 + 26flamq) + 0(a ) (7.9) 

Zs (TrjlPM}) 

is an identity recommended for numerical evaluation. 
7.3 Topological susceptibility 

Using parity-odd density chains, the topological susceptibility Xt in QCD can be 
defined in a manifestly ultraviolet-finite and therefore universally valid way [10]. On 
the lattice there exist different definitions of this type, all of which are expected to 
coincide in the continuum limit. In particular, one can make use of twisted-mass 
density chains and it is then possible, as in the case of the chiral Ward identities 
discussed in the previous subsection, to pass from density chains to spectral projec- 
tors. 

Proceeding along these lines, the expression 

Xt = (1 + 2bRam^) ^:^(I¥{75Pm}'IV{75Pm}) + O(a') (7.10) 
is obtained, which, when combined with eq. (7.9), leads to the formula 

V (Tr|75PM75PM}) 

In principle the mass Mr can be set to any value in eqs. (7.9)-(7.11), but since the 
size of the lattice effects depends on Mr, its value should in practice be chosen with 
some care. One evidently requires that oMr <^ 1 and it is certainly wise to avoid 
the threshold region Mr ~ ttir, where the lattice effects tend to be kinematically 
enhanced. Moreover, a definite prescription that fixes Mr in physical units should 
be adopted when scaling to the continuum limit, as otherwise there is no guarantee 
that the calculated renormalized quantities converge with a rate proportional to . 



8. Concluding remarks 

The condensation of the low modes of the Dirac operator seen in numerical lattice 
QCD provides a most direct piece of theoretical evidence for the spontaneous break- 
ing of chiral symmetry in QCD. Explicit violations of chiral symmetry at momenta 



19 



on the order of the inverse lattice spacing have httle influence on the mode condensa- 
tion, because the mode number is a renormalizable quantity and therefore coincides 
with its continuum limit up to terms that vanish proportionally to a power of the 
lattice spacing. 

The dynamical mechanisms that cause the modes to condense are presently not 
known. It is quite clear, however, that the spontaneous breaking of chiral symmetry 
is not a many-quark collective effect. The mode condensation actually appears to be 
largely insensitive to the sea-quark mass and it seems to persist even when passing to 
the quenched theory. There is thus no relevant back-reaction of the sea quarks and 
theoretical studies of the behaviour of a single quark in presence of representative 
gauge fields may therefore allow the breaking of chiral symmetry to be explained. 

While the computation of the chiral condensate is an obvious application of the 
spectral projector technique introduced in this paper, there are other applications as 
well and the technique is also not limited to a particular lattice formulation of QCD. 
Moreover, it may be useful for studies of the theory at non-zero temperature and of 
QCD-like theories, where chiral symmetry may or may not be spontaneously broken. 

We wish to thank Stefan Sint and Peter Weisz for correspondence on the 0(a)- 
improvement of twisted-mass QCD and Filippo Palombi for his help in producing the 
eigenvalue data on which fig. 3 is based. The gauge-field configurations used for the 
numerical studies were provided by the CLS community [32] . All computations were 
performed on PC clusters at CERN and CILEA. We are grateful to these institutions 
for providing the required resources and their technical staff for assistance. 



Appendix A. SU(3|1) chiral perturbation theory 

As explained in sect. 4, the chiral expansion of the spectral density (and thus of the 
mode number) is obtained by calculating the valence-quark condensate Sva^JTivai) in 
partially quenched chiral perturbation theory [14,15]. We here provide some details 
of this computation, assuming the reader is familiar with chiral perturbation theory 
and partial quenching. 

Following a suggestion of Sharpe and Shoresh [15], we do not include a flavour- 
singlet field in the effective chiral theory. In the sea-quark sector, the chiral expan- 
sions generated by the SU(3|1) chiral lagrangian then literally coincide with those 
obtained in the standard SU(2) theory. In particular, the low-energy constants (such 
as F and S) which already occur in the latter are the same. 



20 



A.l Group generators 

The complex Lie superalgebra of SU(3|1) consists of all 4 x 4 supermatrices XajS 
with vanishing supertrace (see ref. [30], for example). We assume the indices a, /? of 
these matrices to be such that a = 1,2 corresponds to the sea quarks, a = 3 to the 
valence quark and o; = 4 to the ghost (or pseudo-fermion) quark associated to the 
valence quark. 

Our conventions for the generators T", a = 1, . . . , 15, of the algebra are 



T" = {T^y, Str{T"} = 0, Str {TT''} = |c/"^ 
where the non-zero elements of the matrix g"'^ are given by 



(A.l) 



— T 



— T 



V 



-1/ 



1-8 

9-14 
} 15 



(A.2) 



More specifically, T^, . . . are assumed to be generators of the SU(3) subgroup 
acting on the sea and valence quarks, while T^, . . . , T'^'^ mix the ghost with the other 
quarks and T^^ is a diagonal matrix with a non-zero ghost-quark component. 

In the following the Einstein summation convention is adopted for SU(3|1) group 
indices and for Lorentz indices. It is also helpful to introduce the tensors 



which satisfy A;"'=/i'='' = 0. 
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bl5\ 



•ab 



(5"' + 5"'')(5'" + 5' 



,615^ 



(A.3) 



A.2 Chiral effective lagrangian 

The SU(3|1) chiral effective theory is a non-linear cr-model in which the basic field 
U{x) takes values in SU(3|1). As usual the lagrangian 

C = C^^'> + + . . . (A.4) 

is given as a series of terms of increasing dimension. The leading-order term, 

C(^) =-lF^Sti{J^J^}-^BF^Sti{MU'< +M^U}, J^ = U^d^U, (A.5) 
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involves the quark mass matrix M, the pion decay constant in the chiral hmit, F, 
and the parameter B, which is related to the quark condensate through S = BF^. 
The mass matrix is taken to be diagonal, 



where m, mvai and TTivai are, respectively, the masses of the sea quarks, the valence 
quark and the ghost quark. In order to properly quench the valence quark, jfivai will 
later be set to mvai- 

At next-to-leading order, the effective lagrangian reads 

C^^^ = -Lo Str {J^ J, J^JJ - (Li - iLo)Str{J^J^}Str{J,J,} 

- (L2 - Lo)Str{J^J^} Str {J^J^} - (L3 + 2Lo)Sti {J^J^J^J^} 

- 2BLiStv{J^J^}Stv{MU'' + M^C/} - 25^5 Str{ J^([/tM + M^C/)} 

- AB'^LeStiiU^'M + M^U}Sti{U^M + M^U} 

- AB^LrStiiM^U - MU^Stv{M^U - MU^ 

~4B^LsSti{MU''MU'' +M'iUM'iU}-4B^H2StT{M'iM}. (A.7) 

The additional low-energy constants at this order are thus Lq,...,L^ and H2. Note 

that in these expressions we have omitted all terms that do not contribute to the 
valence-quark condensate (such as those related to current correlation functions, for 
example) [16]. 

A. 3 Perturbation expansion 
The chiral expansion of 



M = diag {m, m, mvai , myai}, 



(A.6) 



Sval("^val) 



»nval = '"val ' 



(A.8) 



= -SReC/33 + -.., 



(A.9) 



<7val — 



Smval 



is obtained by substituting 



U = exp {2z^/F} , 



'a 



(A.IO) 
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in the functional integral and expanding all entries in powers of (f). Since the expec- 
tation value in eq. (A. 8) is to be computed at mvai = ^vab one needs to work out 

the Feynman rules only for this case. 

To second order in (j), the leading-order lagrangian reads 

£(2) = \g-\d^rd^,ct>' + M!r<i>'} + U^l - Mi)k'^'cf>'^(f>\ (A.ii) 

where = 2Bm, = 2Bm^^i and 

'Ml if a = 1,2, 3, 

^1 = \ hiMl + if a = 4, . . . , 7, 9, . . . , 12, (A.12) 

. if a = 8, 13, 14, 15. 

The propagator of the meson field is thus given by 

(r(x)0^O)) = 5"'Gi(x,m2) + l{Ml - Ml)h'^^G2{x,Ml), (A.13) 

/rl4 ipx 

All other Feynman rules can be derived straightforwardly from the lagrangian and 
the field a^ai- 

Following common practice, we use dimensional regularization for the loop in- 
tegrals and a modified minimal subtraction scheme for the bare couplings in the 
lagrangian jC^^h In particular, in 4 — 2e dimensions we substitute 

Le = ^^-^!.--+j-lnATr-l + k] (A.15) 



(327r)2 t e 

for the coupling Lg, where 7 = 0.577 . . . denotes Euler's constant, Iq the renormalized 
coupling and p, the renormalization scale. 

A. 4 Finite-volume correction 
The chiral expansion 

SvaK'TT'val) " ^^val ("ival) |y=oo = 

E 



2F2 



{51 (M^,) - 4ffi (IM^ + iMl) + {Ml - Ml) 92 {M^) } + ... (A.16) 
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starts at next-toleading order and involves the momentum sums 

9n{M^) = ^Ej^^J^ - Gn{0,M^). (A.17) 

These are easily calculated numerically when written in the form of rapidly converg- 
ing series of Bessel functions [17]. 



Appendix B. Estimation of the approximation error A 



The computational strategy outlined in sect. 5 assumes that the parameters n, e, M 
and are such that the approximation error A [eq. (5.8)] can be safely neglected. 
In this appendix, we now show how this condition can be met in practice. 

B.l Spectral integral 

Our starting point is the spectral integral representation 

POO 

A= / du {9{M - u) - h{x^)^} v'(ui,m^ (B.l) 

^0 



in which 



d 2M^ 



Note that u'{u>,mq) coincides with the average spectral density of the square root 
of DrrjDm up to a factor V. For illustration, the two functions in the curly bracket 
are plotted in fig. 5 for a typical choice of the parameters. 

In the following, we distinguish three ranges of lo, separated by the limits 

of the transition region around lo = M* (see fig. 5). The parts of the spectral integral 
corresponding to the integration ranges [0,a;_], [ix'_,u;_|_] and [ti;+,oo] are denoted 
by A_, Aq and A+, respectively. 
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60 90 120 

co[MeV] 



180 



Fig. 5. Approximate spectral step function h{x^Y fo'^ = 32, e = 0.01 and = 94 
MeV. The exact step function 6{M — uj) is also shown (grey line), where M and M* 
are related through eq. (B.7). 

B.2 Bounds on A_|_ and A_ 
Noting 



x„ 



±-\/e at 



UJ 



(B.4) 



and recalling the approximation property (5.5) of the minmax polynomial P{y), the 
function in the curly bracket in eq. (B.l) is easily bounded in the case of the integrals 
A±. Since the total number of eigenmodes of Dm Dm is 12V/a'^ and since there are 
at most u{M,mq) eigenmodes with eigenvalues uj'^ < w^, it is then straightforward 
to establish the bounds 



3 V6^ 

< ^• 

- 4 a4 ' 



(B.5) 



|A_| < iy{M,mq) {25 + 0(5^)} 



(B.6) 



These parts of the total error A are thus controlled by the precision 5 of the poly- 
nomial approximation to the step function. 

In ref. [19] it was shown that 6 is an exponentially decreasing function of n^/e. 
The precision can therefore be set to the desired level by adjusting the degree n of 



25 



the minmax polynomial. If e = 0.01, for example, and if a lattice of size 128 x 64^ 
is considered, a sensible choice is n = 32 and eqs. (B.5),(B.6) then imply that 
|A+| < IQ-s and |A_| < 10"^ x z/(M,mq). 

B. 3 Estimation of Aq and the relation of M to M* 

The remaining error component, Aq, is more difficult to estimate than A+ and A_. 
An important point to note is that the density u'{u>,mq) tends to be practically 

constant in the transition region (cf. sect. 6). Most of the error can therefore be 
cancelled by choosing the relation between M and M* to be such that Aq vanishes 
for a constant density. This condition amounts to setting 

and the residual value of the error, 

Ao = f ^ du {e{M -Lo)- h{x^y} {u'{u, mq) - u'{M, m^)} , (B.8) 

•Jul — 

is then of order e. 

An estimation of Ao however requires some information on the a;-depcndcnce of 
the density z^'(a;,mq) in the transition region. For a determination of the expected 
order of magnitude of Aq, chiral perturbation theory may be used at this point and 
a rough bound on the slope of i^'(a;, mq) (and thus on Aq) can also be obtained a 
posteriori through a fit of simulation results for the mode number. Whether Aq is in 
fact negligible with respect to the statistical errors can ultimatly always be checked 
by varying e at fixed 6. 



Appendix C. Lattice parameters and simulation results 

C.l Lattice parameters 

The numerical studies reported in sect. 6 are based on representative ensembles of 
gauge- field configurations for the two- flavour 0(a)-improved Wilson theory (cf. sub- 
sect. 2.2). The ensembles were generated by the authors of ref. [31] and were made 
available to us through the CLS community effort [32]. 
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Table 2. Simulation results for the mode number 



Run 


Lattice 


At 


^cfg 


aM 




^3 


48 X 243 


0.13610 


160 


0.02674 
0.02377 
0.02087 
0.01808 


28.6(4) 
24.0(4) 
19.5(4) 
15.2(3) 


^5 


48 X 243 


0.13625 


160 


0.02549 
0.02234 
0.01923 
0.01616 


26.8(5) 
22.6(4) 
18.5(4) 
14.4(3) 




64 X 32^ 


0.13610 


80 


0.02674 

0.02377 

n nons? 
u.uzuo t 

0.01808 


93.9(12) 
78.2(12) 

48.9(9) 




64 X 32^ 


0.13625 


80 


0.02549 
0.02234 
0.01923 
0.01616 


83.2(10) 
69.6(9) 

56.4(8) 
44.6(7) 




64 X 32^ 


0.13635 


80 


0.02499 
0.02177 
0.01856 
0.01537 


78.7(11) 
66.6(11) 

54.9(9) 
43.8(8) 



In these simulations, the coupling P = Q/g^ was set to 5.3 in all cases and the 
sea-quark hopping parameter /c = (8 + 2mo)~^ to the values quoted in table 2. The 
lattice sizes and the numbers iVcfg of configurations are also given in the table. The 
spacing of the two lattices considered was determined to be 0.0784(10) fm [31] and 
their spatial sizes are thus L = 1.88(2) fm and L = 2.51(3) fm, respectively. 
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C.2 Computation of the mode number 

The mode numbers listed in table 2 were computed stochastically following the lines 
of sect. 5. We used the same minmax polynomial of degree n = 32 in all cases, with 
e set to 0.01, and the number N of pscudo-fermion fields was taken to be 1. With 
these choices, 6 = 4.4 x lO""*, the integral (B.7) evaluates to M/M* = 0.96334 and 
the approximation error A [eq. (5.8)] is estimated to be neglible in our computations. 

The statistical errors quoted in table 2 are in the range from 1 to 2 percent. They 
are practically given by the last term in eq. (5.4), which explains why, with half the 
statistics on the larger lattices, approximately the same relative accuracy is obtained 
on all lattices. Moreover, the errors could be further reduced, by a factor 2 at least, 
by increasing the number N of pseudo-fermion fields. 

C.3 O (a) -improvement and renormalization at (3 = 5.3 

The coefficients of the 0(a) counterterms in the quark action [3] and the improved 
axial current [4] were set to the non-perturbatively determined values Cgw = 1.90952 
[6] and ca = —0.0506 [33], respectively. We computed the renormalized quark mass 
via the PCAC relation, using the improved axial current, but neglected the O(amq) 
corrections in eq. (2.8) since 6^ — 6p = —0.00104(6) x (/q + 0((7q) is very small [35]. 

AlthoTigh a different improvement scheme was adopted in ref. [8], it is possible to 
deduce the one-loop formulae 



from the results published there and in refs. [34,35]. So far is only known in per- 
turbation theory and we thus used the one-loop estimate = —0.626 in eq. (3.11). 
Noting arric = —0.33560(5), the subtracted bare mass amq is smaller than 0.01 at 
all values of n considered. The calculated O(amq) corrections to Mr are therefore 
at most 0.6% and the corrections to the ratio (7.9) will normally be negligible. 

The renormalization factors Za = 0.778(10) [36] and Zp = 0.543(8) [37,38] needed 
to pass from the bare masses m and M to the renormalized masses and Mr in 
the MS scheme at 2 GcV have been computed non-perturbatively. As can be seen 
from table 1, the renormalized sea-quark mass ranges from about 13 to 46 MeV on 
the lattices considered. The values of aM in table 2 have, incidentally, been chosen 
such that Ar = (M^ — m^)^/^ approximately assumes the values 70, 85, 100 and 115 
MeV at all sea-quark masses. 




-i - 0.111(4) X 5^ + 0(5^), 



(C.l) 



bR = -0.031(8) xgl + 0{g^), 



(C.2) 
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